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I.  INTRODUCTION 


We  examine  here  a  method  for  reconstructing  an  x-ray  attenuation  function  from 
measurements  of  its  integrals  in  the  case  where  the  projection  data  is  sparse.  In  practice, 
data  from  as  few  as  5  to  20  projections,  with  perhaps  only  25  measurements  per  projection, 
may  be  all  that  is  available.  This  is  the  case,  for  example,  in  industrial  nondestructive 
testing,  where  the  object  whose  density  cross-section  we  are  attempting  to  determine  is  in 
a  rapid  state  of  flux.  In  this  event,  there  is  very  little  time  available  in  which  to  gather 
data.  Farther,  what  little  data  there  is  available  is  almost  certainly  degraded  by  noise  and 
possibly  blurring  due  to  the  motion  of  the  object. 

One  should  note  that  the  standard  reconstruction  algorithms  currently  in  use  in  most 
medical  CAT  scanners  utilise  180  views  (or  projection  angles)  and  yield  very  poor  resolution 

when  the  number  of  views  is  substantially  reduced*.  One  approach  which  has  proven  useful 
when  the  available  projection  data  is  sparse  is  the  maximum  entropy  algorithm  MENT  of 

A 

Vinerbo.  In  this  method,  one  attempts  to  compute  the  attenuation  function  having  maximum 
entropy  (intuitively,  one  might  say,  the  function  yielding  the  picture  with  the  least 
2  3 

information  content  ’  )  and  which  is  also  consistent  with  the  measurements.  One  of  the 
difficulties  with  this  approach  is  that  Vinerbo  looks  for  a  solution  in  the  class  of  functions 

continuous  in  the  scanning  region,  while  it  has  been  shown  by  Klaus  and  Smith4  that  the 
solution  of  this  problem  must  be  piecewise  constant  throughout  the  scanning  region.  Also,  the 
MENT  algorithm  as  now  implemented  does  not  take  advantage  of  a  priori  information  known 
about  the  object  being  scanned.  We  address  both  of  these  problems  in  the  current  work. 
The  most  significant  difference,  however  between  MENT  and  the  current  work  is  that 
MENT  requires  that  the  density  of  the  object  being  scanned  tends  to  sero  as  the  edges  of 
the  scanning  region  are  approached.  If  this  requirement  is  not  met,  MENT  may  produce  a 
rmnUuly  unrecognisable  image-  Naturally,  the  only  way  to  guarantee  that  such  an  "edge 
condition"  is  satisfied  in  practice,  is  to  require  minimum  x-ray  source-to-phantom  and 
phantom- to- detector  distances,  where  the  density  of  the  surrounding  medium  is  known  a 
priori.  This  condition,  then  imposes  a  restriction  on  the  maximum  site  of  an  object  to  be 
scanned  by  a  given  device,  a  sise  smaller  than  that  which  the  machinery  itself  would 
otherwise  allow.  The  current  method  removes  this  restriction  entirely.  We  have  included  an 
example  of  a  reconstructed  density  profile  with  nonsero  density  near  the  edge  of  the 
scanning  region,  obtained  by  the  new  method.  As  a  comparison,  we  include  the  MENT 
reconstruction  obtained  using  the  same  data. 


*C.K.  Zoltani,  K.J.  White  and  R.P.  Kruger,  "Result  of  Feasibility  Study  on  Computer  Assisted 
Tomography  for  Ballistic  Applications",  ARBRL-TR-02513,  Aberdeen  Proving  Ground. 
Maryland,  ADA  133  214,  August,  1983. 

2 

G.  Vinerbo,  "MENT:  A  Maximum  Entropy  Algorithm  for  Reconstructing  a  Source  from 
Projection  Data",  Comp.  Graph,  and  Image  Proc.,  Vol.  10,  pp  48-68,  1979 

*RD.  Levine  and  M  Tribus,  The  Maximum  Entropy  Formalism.  MIT  Pi  ess,  Cambridge.  1978 

4M  Klaus  and  R.T  Smith,  "A  Hilbert  Space  Approach  to  Maximum  Entropy  Reconst  ruction". 
Math  Meth  in  the  Appl  Set.,  to  appear 


In  [4J,  Klaus  and  Smith  demonstrated  that  by  subtracting  the  entropy  from  a 

penalty  term  consisting  of  the  residual  error  in  matching  the  measurements  and  then 

minimising  this  augmented  cost  functional  over  a  weakly  compact  set,  one  obtains  a 

2  2 
well-posed  optimisation  problem  in  L  +(D),  where  D  is  the  scanning  region  and  l/+(D)  is 

the  space  of  non -negative  square  integrable  functions  on  D.  That  is,  there  exists  a  unique 

solution  of  the  problem  and  the  solution  depends  continuously  (weakly  in  L2+(D)  )  on  the 

measured  data.  It  was  also  shown  that  the  solution  of  this  problem,  for  the  special  case  of 
the  weakly  compact  set  which  we  consider  here,  must  be  piecewise  constant.  Our  problem  is 

A 

set  in  the  larger  function  class  L  (D),  rather  than  in  the  class  of  piecewise  constant 
functions,  however,  since  the  geometry  of  the  sets  on  which  the  solution  is  constant  is  so 
complicated  as  to  make  trying  to  find  the  exact  solution  of  the  optimisation  problem 
prohibitive. 


In  the  present  work,  a  solution  of  the  optimisation  problem  is  sought  over  a  finite 
o 

dimensional  subspace  of  L  (D).  Thus,  the  infinite-dimensional  optimisation  problem  is  reduced 


to  one  over  Rn.  A  simple  discretisation  of  D  is  made  and  a  space  of  finite  elements  is 
defined  on  this  grid.  We  restrict  ourselves  here  to  basis  functions  which  are  piecewise 
constant  or  products  of  piecewise  linear  functions  in  each  variable.  It  is  shown  here  that  as 
the  finite  element  mesh  sise  tends  to  sero,  the  sequence  of  solutions  of  the  finite 

A 

dimensional  optimisation  problems  will  converge,  in  L  (D)  to  the  solution  of  the  original, 
infinite  dimensional  problem. 


The  finite  dimensional  optimisation  problems  are  somewhat  large  in  practice  (usually 
at  least  300  to  1,000  variables,  for  minimal  resolution)  but  they  are  relatively  simple  in 
form.  Although  the  cost  functional  is  quite  nonlinear,  the  constraints  are  very  simple  bound 
constraints. 

For  our  numerical  experiments,  we  have  used  the  MINOS  nonlinear  optimisation 

package  of  Stanford’s  Systems  Optimisation  Laboratory6’  .  The  current  version  of  MINOS 
uses  a  reduced  gradient  scheme  together  with  a  quasi-Newton  method. 


II.  THE  OPTIMIZATION  PROBLEM 


We  consider  the  problem  of  reconstructing  an  x-ray  attenuation  function  in  two 
dimensions  from  measurements  of  its  integrals.  The  x-rays  are  assumed  to  be  collimated  so 
that  a  parallel  beam  scanning  geometry  is  obtained.  Assume  that  J  x-ray  sources 


6P.  Gill,  W.  Murray  and  M  Wright,  Practical  Optimisation.  Academic  Press.  London,  1981 

®B.A.  Murtaugh  and  MA.  Saunders,  "Large-Scale  Linearly  Constrained  Optimization".  Math 
Prog.,  Vol.  14,  pp.41-72,  1978 


are  arranged  on  a  circle  centered  at  the  origin  and  containing  the  region  to  be  scanned,  D  c 

R2.  Suppose  that  the  center  lines  of  the  x-ray  projections  make  angles  ay  j  —  1,...,J, 

with  the  positive  x-axis.  We  assume  that  we  have  available  measurements  of  the  x-ray 
attenuation  along  strips  A^,  m  ~  l,...,M(j),  which  are  parallel  to  the  center  line  of  the 

x-ray  projection  from  the  j1^  source,  for  each  j  —  1,...,J.  The  strips  are  such  that  for  each 

j  "  1  r”4. 

M(  j) 

Dc  u  A.  . 

jm 

m—  1 


o 

The  entropy  of  f,  for  f  €  L  (D),  is  defined  to  be 
*(f)  -  -  |  f(x,y)  |  In  |  |  f(x,y)  |  A  )  dx  dy  , 


where  A  is  the  area  of  D. 

Let  f q  be  the  actual  x-ray  attenuation  function.  We  will  assume  that  there  is 
measured  projection  data  available  in  the  form 


<>  -  Vf0):- 


fQ(x,y)  dx  dy  , 

*  « 


m  -  1,  -Mi) ;  j  -  1.  4 


We  add  to  the  entropy  a  penalty  term  corresponding  to  the  residual  error  in  meeting 
the  measurements  (2.2).  Define  the  functional 


G(f)  -  - 


n(f)  +  T  £  (G^  -  s^f)  ]* 


where  7  >  0  is  the  penalty  parameter. 

The  object  then,  is  to  minimize  (2.3),  subject  to  some  appropriate  constraints  on  f. 
This  should  yield  the  picture  with  the  least  information  content  which  matches  the  measured 
data  to  within  a  specified  error.  (Note  that  the  residual  error  is  reduced  by  taking 


7 

larger  values  of  the  penalty  parameter  'i  .]  We  then  consider  the  optimization  problem 

inf  G(f)  ,  (2.4) 

f  6  E 

where  E  is  a  weakly  (sequentially)  compact  subset  of  L^+(D).  By  weakly  (sequentially) 

compact,  we  mean  that  every  sequence  in  E  has  a  weakly  convergent  subsequence  whose 
weak  limit  lies  in  E.  Of  course,  (xQ)  weakly  convergent  to  x  means  that  the  sequence  of 

A 

real  numbers  {(xn,y)}  converges  to  (x,y)  for  every  y  in  L  (D),  where  (’,  )  denotes  the  usual 

L  inner  product.In  practice,  E  will  contain  a  priori  information  known  about  the  object 
being  scanned  (e.g.  E  may  specify  maximum  and  minimum  densities,  information  on  the 
support  of  the  attenuation  function,  etc.). 

In  [4],  it  was  shown  that  the  optimization  problem  (2.4)  is  well-posed.  That  is,  there 
exits  a  unique  solution  for  any  given  set  of  measured  data  and  the  solution  depends 

Q 

continuously  (weakly  in  L  (D)  )  on  the  measured  data.  The  key  to  the  proof  of 

well-posedness  hangs  on  the  following  two  results,  which  will  also  be  needed  to 
demonstrate  the  convergence  of  our  numerical  method. 

Lemma  1.: 

o 

The  cost  functional  G  is  continuous  on  L  (D). 

Lemma  £ 

2 

G  is  a  strictly  convex,  weakly  lower  semicontinuous  functional  on  L  +(D).By  weakly 

2 

lower  semicontinuous,  we  mean  that  if  (fn)  converges  weakly  to  f  6  L  +(D),  then 

G(f)  *  J_im  G(fn). 

n  — ►  oo 

For  the  proof  of  Lemmas  1  and  2,  the  interested  reader  is  referred  to  (4]. 

It  is  also  shown  in  (4),  that  the  solution  of  (2.4)  for  the  special  case  of  E 
considered  here  is  constant  on  each  cell  of  D,  where  by  a  cell,  one  means  the  interior  of  a 
maximal  region  which  is  not  crossed  by  any  of  the  strip  boundaries.  This  characterization  of 
the  solution,  while  of  theoretical  interest,  is  of  little  help  in  actually  solving  the 
optimization  problem,  because  of  the  very  complicated  geometry  involved.  Rather  than  using 

this  characterization  of  the  solution,  then,  we  consider  a  method  of  solution  utilizing  the 
o 

L  (D)  structure. 

7R  Fletcher,  Practical  Methods  o£  Optimization.  Vol.  2,  John  Wiley,  New  York,  1980. 
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One  should  note  that,  due  to  the  severe  nonlinearity  of  the  object  function  G(f),  it 
is  rather  difficult  to  approach  (2.4)  as  an  infinite  dimensional  problem  using  the  calculus  of 

O 

variations.  We  therefore  approach  the  problem  via  the  Ritx  method  .  For  each  h  >  0,  let 

L  A 

S  denote  a  finite  dimensional  subspace  of  L  (D),  satisfying  the  following  approximation 
property  GSven  e  >  0,  one  can  choose  h  sufficiently  small  so  that  the  best 

approximation  f*1  to  f  out  of  (Le.  the  projection  of  f  onto  S*1)  satisfies 


Consider  now  solving  the  optimisation  problem  (2.4)  over  S^,  that  is, 

inf  G(f)  ,  (2.6) 

f  e  Ensh 

for  a  sequence  of  values  of  h,  tending  to  sero.  We  can  now  prove  the  following  convergence 
result. 

Theorem: 

Let  {  S*1  |  h  >  0  }  be  a  family  of  finite  dimensional  subspaces  of  L2(D)  satisfying 
the  above  approximation  property.  Then  for 


#•  h*-G<0-  inf  ,  <H0. 

f  €  EnSh 

*  *  *  2  * 

/ijj  converges  to  /i  and  fjj  converges  weakly  (  in  L  (D)  )  to  f  ,  as  h  tends  to  0, 

where 


£rsof: 


! i *  -  G(f*)  -  inf  G(f). 

f  e  E 


2 

Since  G  is  weakly  lower  semicontinuous  and  strictly  convex  on  L  +(D)  and  E  is 
weakly  compact,  there  must  exist  a  unique  minimizer,  of  G  over  E  fl  S*\  for  each  h  > 

0.  Let  f*1  denote  the  best  approximation  to  f  out  of  S*1.  Since  G  is  continuous  on  L2(D) 

(Lemma  1)  and  the  family  of  subspaces  S^1  satisfies  the  above  approximation  property,  one 
may,  given  <  >  0,  choose  h  >  0  sufficiently  small  so  that 


g 

I.M  Gelfand  and  S.V.  Fomin,  Calculus  2I  Variations.  Prentice  Hall,  Englewood  Cliffs,  New 

Jersey,  1963. 


(2.7) 


G(fh)  <  G(f*)  +  *. 

Note  that 

ti*  -  G(f*)  -  inf  G(f)  £  inf  G(f)  -  G(fh*).  (2.8) 

f  €  E  f  €  EnSh 

However,  since  f*1  €  S*\ 

G(fh*)  -  inf  G(f)  £  G(fh). 
f  e  Ensh 

From  (2.7)  and  (2.8),  it  then  follows  that 

H*  Z  G(fh*)  £  G(fh)  <  G(f*)  +  t  -  m*  +  <> 

*  * 

and  hence,  G(f^  )  converges  to  n  ,  as  h  tends  to  0.  Note  that  this  alone  is  not  sufficient 

to  say  that  f^  converges  to  f  in  L  +(D).  Since  one  is  interested  in  approximating  the 

attenuation  function  and  not  just  its  associated  entropy,  this  last  question  is  crucial.  This 
can  be  resolved  as  follows. 

Let  (hn)  be  a  sequence  of  positive  real  numbers  tending  to  0  and  let 


-  G(f  ) 


inf  h  G(f). 
f  €  EPS  n 


Then  {  f„  }  is  a  minimizing  sequence  for  G(f)  over  E,  from  the  preceding.  Since  E  is 

*  * 
weakly  (sequentially)  compact,  {  f„  }  must  have  a  weakly  convergent  subsequence  {  f  } 

** 

whose  weak  limit  f  lies  in  E.  It  follows  by  the  weak  lower  semicontinuity  of  G  that 


G(f**)  ^  ilm  G(fn  *) 

k  — *  oo  k 


G(f  ). 


**  * 

whence  f  —  f  ,  since  G  has  a  unique  minimizer  in  E.  Lastly,  since  E  is  weakly  compact, 
*  * 

and  {  fn  }  C  E,  every  subsequence  of  {  }  has,  in  turn,  a  subsequence  converging 

*  *  * 
weakly  to  f  .  Thus,  {  f  }  must  itself  converge  weakly  to  f  . 


12 


One  can  see  from  the  theorem  that,  by  solving  the  optimization  problem  (2.6)  for  a 

finite  dimensional  subspace  Sa  of  L  (D),  we  do  indeed  obtain  an  approximation  to  the 
solution  of  the  infinite  dimensional  problem  (2.4).  In  the  next  section,  we  give  examples  of 
two  finite  element  spaces  in  which  we  have  solved  (2.6).  We  have  included  examples  of 
images  reconstructed  with  this  method  using  extremely  sparse  data. 

m.  FINITE  ELEMENT  APPROXIMATION. 

In  this  section,  we  show  how  to  implement  the  method  described  in  the  last  section 

o 

by  using  a  space  of  finite  elements  as  the  approximating  subspace  of  L  (D).  We  give  several 
examples  of  simple  finite  element  spaces  and  reconstruct  several  images  using  one  of  these 
spaces. 

First,  superimpose  a  fixed  rectangular  mesh  on  D  ■  (—1,1)  x  [-1,1],  with  uniform 
mesh  size,  h  —  1/n  in  both  the  x  and  y  directions,  for  simplicity.  Consider  using  products 

of  piecewise  linear  functions  in  x  and  y  as  the  finite  element  space  S^1.  A  basis  for  S*1  is 
given  by 

^(x.y)  “  ^(y)  ,  k  -  l,...,(2n+l)2, 

where 

1  -  [(k-l)-(k-l)  mod  (2n+l)]/(2n+l)  -  n, 
i  -  k  -  (l+n)(2n+l)  -  n  -  1 

and 

o  ,  t  £  ( j -l ) h ,  t  £  ( j  +  l)h 
<f> j(t)  -  [t-(j-l)h]/h  ,  (j-l)h  £  t  £  jh 

[ ( j+1 )h-t ] /h  ,  jh  £  t  £  ( j+1 )h . 

It  is  very  reasonable  to  expect  that  in  practice,  one  should  know  a  priori  the 
minimum  and  maximum  densities  of  the  object  being  examined.  Therfore,  we  have  defined  a 
simple  constraint  set  by 

E  —  (f€  L“(D)  |o<a^fiab<oo,  a.e.  and  f  =  0  a.e.  in  F2  s  D) 

The  attenuation  function  f  is  then  approximated  in  S*1  by 

N 

f(x,y)  =  I  c,  <Mx,y), 
k=  1  K  K 

2 

where  N  —  (2n  +  l)  and  the  coefficients  c^  are  chosen  as  the  solution  of  the  finite 
dimensional  optimization  problem 


(3.1) 


inf  G  (  2  c.  ^.(x,y)  ) 


c  6  R 


,N  k-1 


subject  to 


0  <  a  <  I  c.  ^.(x,y)  £  b. 


With  the  simple  choice  of  the  finite  element  space  Sn  and  due  to  the  form  of  G, 
(3.1)  -  (3.2)  reduces  to 


inf  {  -V  (  ^  ck  ^k(x,y)  )  +  'f  ^  I  Gi 

o  6  rn  k-‘  i  ’ 


-  ck  Sjm<  *k<x^  )  ]2  } 


subject  to 


5  c^  ^  b,  k  — 1,2,...,N. 


The  finite  dimensional  optimization  problem  (3.3)  -  (3.4)  can  be  solved  numerically, 
for  sufficiently  small  values  of  N,  using  existing  optimization  software  packages.  [For  our 

_ i  -  ~  i_ _  _ ■  ,l.  _ l _ 6  i  m...  .  u  •  _ l.  J _ _ l—  :  r _ i 


examples,  we  have  used  the  MINOS  package  .]  Note  that  this  can  be  done  only  if  one  has  a 
simple  method  of  computing  values  of  the  object  function  for  each  updated  value  of  Cj,...,Cj^. 

Because  of  the  form  of  (3.3),  one  needs  to  compute  the  values  of  Sjm(  ^(x,y)  )  for  each 

j,  m  and  k  only  once  for  given  projection  and  discretization  geometries.  These  values  can 
then  be  stored  in  a  (somewhat  large)  disk  file  and  simply  read  back  each  time  new  values 
of  Cj,...,Cj^  are  computed  (i.e.  in  each  pass  of  the  optimization  scheme).  We  have  devised  an 

algorithm  which  can  compute  the  integrals  Sjm(  0j.(x,y)  )  exactly  for  any  given  projection 

and  discretization  geometries.  It  is  then  a  simple  matter  to  compute  the  residual  error 


1  Gjm  '  kti  Ck  V  )  J 


at  each  pass  of  the  optimization  scheme. 


Note  that  due  to  the  nonlinearity  of  the  entropy  term,  we  need  to  recompute 


n(  I  C,  <Mx,y)  ) 
k-1  K  K 


are 


for  each  new  set  of  values  Although  the  geometry  is  simple,  the  integrands 

rather  complicated.  One  can  write  the  entropy  as  an  iterated  integral,  the  first  of  which  is 
evaluated  explicitly.  The  second  integration  is  then  performed  numerically,  each  time  a  new 
set  of  values  of  Cj,...c^  is  computed.  In  the  examples,  we  have  used  an  adaptive  quadrature 

routine  for  this. 

The  question  of  how  to  determine  an  optimal  value  of  the  penalty  parameter  7  is 
as  yet  unresolved.  In  practice,  7  must  be  sufficiently  large  so  that  the  residual  error  in 
meeting  the  measurements  (3.5)  is  small.  At  the  same  time,  7  must  not  be  so  large  that  the 
penalty  term  completely  swamps  out  the  entropy.  Computational  experience  has  been  our 
only  guide  thus  far,  for  determining  an  acceptable  value  of  7.  More  work  remains  to  be 
done  in  this  direction. 

As  yet,  nothing  specific  has  been  said  about  how  to  choose  the  mesh  size  h.  In  the 
theorem  of  the  last  section,  we  have  seen  that  as  h  tends  to  0,  the  solution  of  (3.3)  - 
(3.4)  tends  to  that  of  the  infinite  dimensional  problem  (2.4).  In  practice,  however,  one  must 
choose  a  value  of  h  which  will  yield  a  reasonable  resolution  and  yet  for  which  the  computer 
time  is  not  excessive.  There  is  of  course,  also  an  important  restriction  placed  on  h  by  our 
desire  to  be  able  to  choose  7  sufficiently  large  so  that  the  residual  error  (3.5)  is  within 
tolerances.  In  order  that  one  might  make  (3.5)  arbitrarily  small,  it  is  necessary  that  the 
mesh  be  sufficiently  fine  that  the  system  of  equations 

N 

Gjm  ”  JSj  ck  V  *k(x>y)  )  -  m  -  ^  «  ;  j  -  i.-.J 

is  consistent  (i.e.  there  must  be  at  least  one  solution).  For  the  examples  which  we  have 
considered,  we  have  screened  prospective  values  of  the  mesh  size  by  numerically  checking 
the  system  of  measurement  equations  for  consistency.  We  then  increased  the  value  of  the 
penalty  parameter  until  a  minimizer  with  an  acceptable  residual  error  was  found. 

As  an  illustration  of  the  reconstruction  algorithm  presented  above,  we  consider  a 
phantom  consisting  of  three  identical  cylinders,  all  of  density  1,  placed  at  arbitrary 
positions  within  the  scanning  region,  where  the  background  density  in  the  scanning  region  is 
.1  (see  Figure  1  for  a  cross-sectional  view).  Note  that  one  of  the  cylinders  is  rather  close 
to  the  edge  of  the  target  grid.  For  the  reconstruction,  the  phantom  is  set  within  a  square 
grid  30  unit  cells  on  a  side. 


g 

The  computer  code  GOLEM  was  used  to  produce  the  x-ray  absorption  data  from 
the  phantom,  mathematically.  Data  was  produced  from  25  detectors  for  each  of  5  projection 
angles,  arranged  uniformly  ar  ound  the  target.  The  reconstruction  grid  spans  the  width  of 
the  25  detectors,  as  seen  in  Figure  1.  The  31x31  grid  is  used  as  the  finite  element  grid 
and  in  Figures  2  and  3  can  be  seen  inside  a  51x51  grid,  where  points  outside  the  31x31 
grid  are  given  density  zero. 


9 

MD.  Altschuler,  T.  Chang  and  A.  Chu,  "Rapid  Computer  Generation  of  3-D  Phantoms  and 
Their  Cone-Beam  X-ray  Projections",  Medical  Image  Processing  Group  Technical  Report 
MIPG16,  State  University  of  New  York  at  Buffalo,  Nov.  1978. 


SAMPLE  PROBLEM 


XC  YC  RflO  DENSITY 


0ACK8ROUNO 

0.0  0.0 

15.0 

0.1 

1 

-7.0  -7.0 

3.6 

1.0 

2 

0.0  10.7 

3.6 

1.0 

3 

7.2  0.0 

3.6 

1.0 

Figure  1.  Cross 

-Sectional  View  of 

Phantom 

In  Figure  2,  we  give  the  reconstructed  density  plot,  using  our  new  method,  for  the 
phantom  in  Figure  1.  Notice  that  the  three  cylinders  stand  out  clearly  from  the  background. 
In  Figure  3,  we  give  the  MENT  reconstruction  using  the  same  data  as  that  used  for  Figure 
2.  Of  course,  this  is  completely  unrecognizable,  due  to  the  large  mass  near  the  edge  of  the 
target  grid. 

A  more  extensive  collection  of  examples  is  given  in  [10].  More  of  the  technical 
details  of  the  reconstruction  are  given  there,  along  with  a  more  detailed  comparison  with 
MENT. 
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Figure  2.Reconstructed  Density  Plot  using  New  Method 


Figure  3.  Reconstructed  Density  Plot  Using  ME  NT 
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